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We present the renormalized Taylor expansion of the chiral condensate with chemical potential. 
The relation with the Taylor coefficients of meson susceptibilities and the quark mass dependence 
of the quark number susceptibilities are obtained through chiral Ward identities and Maxwell re- 
lations. The continuum limit is obtained in quenched and two flavour QCD. In Nf = 2 QCD, the 
quadratic response coefficients (QRCs) of the chiral condensate to chemical potential indicate that 
pion fluctuations decrease at finite chemical potential toward the critical end point. We discuss how 
these observables can be used to test other features of the QCD phase diagram. 
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I. INTRODUCTION 



Studies of hot dense matter formed in heavy-ion colliders such as the RHIC [l| need to be interpreted using 
theoretical predictions obtained from lattice computations. Most of these are performed at zero baryon density. 
However, it seems likely that the actual experimental conditions involve small but non-vanishing baryon chemical 
potential. In order to deal with this one has to consider the stability of various lattice predictions to the presence of 
£N) ■ a small chemical potential, fj,. 

This was once considered an insoluble problem in lattice QCD. The reason is that direct lattice Monte Carlo 
simulations at finite /x are impossible since the weights in the partition function are not positive definite. However, 
developments in recent years 0, HI 0> mean that QCD at finiteu is no longer out of reach on the lattice. In this 
paper we investigate Taylor expansions in the chemical potential |j, |fj, |j| , which have recently been used to compute 
the equation of state at finite \x 0,0. 

The first investigations of correlators at finite chemical potential was performed in ■ Our approach differs from 
theirs in that we do a direct Taylor expansion of operators, rather than of parameters in fitting functions. There is a 
O ' conceptual and practical gain. Putting chemical potentials on quarks can break the flavourxCP symmetry of QCD 
in many different ways, possibly affecting every observable in distinct fashion. However it turns out that there are 
^-h enough relations between operators and their derivatives that, among all the observables related to chiral dynamics 
in QCD, there is only one set of independent Taylor coefficients. Taylor expansions of the operators utilize these 
relations with generality and economy. 

There turn out to be three independent parameters of interest in the chiral sector — a single linear response 



coefficient of the chiral condensate to \i, and two quadratic response coefficients (QRCs). The linear coefficient 
vanishes by symmetry in an expansion around fi = 0. Thus, the two chiral QRCs encapsulate the physics arising from 
the influence of baryon dynamics due to fj, ^ on chiral fluctuations. Since the phase diagram of QCD at small /j, is 
built up through the interplay of chiral and baryon dynamics 0, E3> EI i the QRCs provide direct information on 
this phase diagram. 

Here we present the first systematic study of the renormalized Taylor expansion of the chiral condensate in QCD 
up to second order in the chemical potentials jl2j . We define renormalization schemes which are simple to implement, 
and in which the Taylor expansion can be extrapolated to the continuum. We find the continuum limits at both high 
and low temperatures in quenched and dynamical Nf — 2 QCD. 

In section 2 we write out generic Taylor expansions in the chemical potential, before specializing to pure glue 
operators, the chiral condensate and meson susceptibilities. The renormalization schemes for the condensate, Maxwell 
relations and chiral Ward identities are also dealt with in this section. The actual lattice simulations and results are 
in section 3, including the complete specification of the renormalization schemes and results above and below T c in 
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the quenched theory. This section also contains our results for QCD with dynamical quarks. Discussion of the results 
and the constraints that the QRCs can place on the QCD phase diagram can be found in the concluding section 4. 

II. TAYLOR EXPANSIONS 

The partition function of QCD at temperature T and chemical potentials fif for each flavour / of quarks is 

Z({m f }, T, { M/ }) = f VUc- s ^ J] Det M(m f , T, fif), (1) 
3 f 

where m/ is the mass of quark flavour /. S(T) denotes the gluon part of the action and M is the (lattice discretised) 
Dirac operator. The expectation value of an operator J is defined as 

(J) ee J({m f }, T, {fi f }) = 4 / VUJe- 3 ^ [] Det M(m f , T, fi f ), (2) 

3 f 

where J may or may not contain explicit dependence on the parameters. Direct simulations are possible only along 
certain subspaces of the full parameter space, such as Re/i/ = or \x u = —/id- However, Taylor expansions around 
fif = have proved to be a good method of studying the physics for general fif. In the remainder of this paper we shall 
work mainly with the two flavour case / £ {u,d} and with degenerate quark masses m u = rrid = m. We shall freely 
go from the parameter space expressed in terms of fi u and fid to the isoscalar chemical potential fiQ = (fi u + fid)/ 2 
and the isovector fi 3 = (fi u — (id)/2- These are simply related to the baryon and electric charge chemical potentials. 
Since we perform Taylor expansions around the point fif = 0, all angular brackets in the rest of this paper shall mean 
expectation values at vanishing chemical potential. These are computable by the standard methods of lattice gauge 
theory. 

We recall a few points about the expansion of the partition function 0, Q . The fi / appear only in the Fermion 
determinant. If a matrix M(x) depends on a parameter x, then the identity T)etM(x) = exp[tr lnM(x)] yields, 
[DetM(x)]' = Det M(x)tr M _1 M', where primes denote derivative with respect to the parameter x. Furthermore, 
the identity M _1 M = 1 yields [Af -1 ]' = -M^ X M'M~ X . This is enough to determine the Taylor coefficients of the 
partition function. The first order coefficients are proportional to the quark number densities — 



dZ 
dfif 



= ZZ 1: where Z x = (d> = (tr M^M}), (3) 

M/=0 



which is flavour independent since m u — md- When we need to keep track of the flavour we shall use the notation 
0\ and obvious generalizations for higher derivatives. The second order coefficients are 



d 2 Z 



= ZZ^ f9) = Z{O u + 2 Sf g ) where (0 2 ) = -{tr Mj 1 M'fMj 1 M' f - trM^M/"). (4) 

n f =o 



The notation used above is On = 0\0\. Flavour symmetry gives two independent components — = Z^ d = 
(0 2 + On) and Z^ 4 = (On). Expansions for isoscalar = fi d = fj, ) and isovector (ji u = —fid = ^3) chemical 
potentials are — 



71 rj, s 7i „ n nw / I + (20n +O 2 )^ + 0( M 4) 

Z(m, T, M „, W ) = Z(rn, T, 0, 0) X | x + + Q 



(5) 



The terms of odd order vanish by CP symmetry. Diagrammatic representations of these operators have been written 
down [l3l | and can be extended to other operators, as we outline later. 

The double series expansion in terms of (j, u and fid of the pressure, P = T(\ogZ)/V, has as coefficients the 
quark number densities, rif = (T/V)(Oi) — 0, and the quark number susceptibilities. The two independent linear 
susceptibilities 0] are the diagonal ones Xuu — Xdd = (T/F)(0 u + 2 ) arising from taking two derivatives with 
respect to the same chemical potential, and the off-diagonal Xud = (T/V)(On). In writing these expressions we have 
used the fact that the expansion is around fif = 0. Higher order susceptibilities have also been used Q|. 

Similar Taylor expansions can be written for any operator J . Here we write this out to second order. Using the 
notation for the first derivative of J with respect to any explicit appearance of/// in it, and J^9) f or the explicit 
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double derivative with respect to /// and n g we can write the first two Taylor coefficients in the expansion- 



J 



(/) _ 



dJ 
dfif 



M/=0 



J. 



(fa) - 



d J 



=0 



The corresponding Taylor expansion is 



J{m,T^ u ,^L d ) = J + j[ u) (i u + j[ d) + 



= ( ( On + S fg 2 )J + 1 [J 



r(f) 



J 



(a) 



J 



(fa) 



(0) 



1 r 
2 



(4 UU) - JoZt } ) A + (4 dd) - JoZ { 2 dd) ) + 2 ( 



(ltd) 



— JqZ 2 - 



(ud) 



(J-u^d 



(7) 



where Jo = J(m,T, 0,0). CP symmetry cannot be used to simplify the expansion in general, since the operator may 
not transform nicely under CP. Special simplifications occur for purely gluon operators, as we discuss later. 

A. Pure gauge operators 

Since pure gauge operators have no explicit dependence on the chemical potentials, 



Ji = (OiJ) , and J 



(fa) 



The corresponding Taylor expansion is 

J(m,T,fi u ,fi d ) = J + J\{nu + (id) + ^ ( 



J 



(uu) 



— JqZ 2 



(uu) 



o 11 + s fg o 2 )J 



0^ + /4) + 2(j 2 ( " 



d) -j 7 {ud) 



(9) 



In general the linear term vanishes exactly only when /i = 0. However, for CP even operators, such as real parts of 
traces of products of links, the odd order terms can be shown to vanish by summing over CP orbits before performing 
the remaining part of the path integral. 



B. Chiral Condensate 



From the flavoured condensates (uu) and (dd), one can build the isoscalar and isovector condensates which are, 
respectively, 



C s (m,T,(i u ,(i d ) = - [(uu) + (dd)] 



and 



1 



C v (m,T,fi u ,fi d ) = - [(uu) - (dd)] 



(10) 



For equal mass quarks, C'v vanishes at vanishing chemical potential. This vector flavour symmetry is not broken by 
/io since a sign flip in /iq carries the quark of each flavour to its own antiquark. As a result, Cy vanishes at all /io, 
and the series expansion for Cs(fJ-o) is even. A sign flip in ^3 carries a quark of a given flavour to the antiquark of 
the opposite flavour, thus mixing the vector flavour symmetry with CP. The same arguments can still be followed 
through to show that Cy vanishes to all orders in the expansion around /X3 = and CsifJs) is even. In the remainder 
of this paper we consider only Cs- 

Recall that the quark condensate can also be written as 



1 f dhgZ dlogZ 



(11) 



m u —7nd—m 



where V4 is the 4- volume of the lattice [l4T | . Spontaneous breaking of chiral symmetry is signaled by a non-zero value 
of the condensate in the limit m = 0. We shall not be concerned with this limit in this paper, and can consequently 
neglect the known subtleties of taking the limits m — > and V4 — > 00. We shall however, check that in actual 
numerical work the volumes are large enough to obtain correct physical results. 

Since the mass is subject to renormalization as the cutoff changes, one must bother about cutoff dependence of 
the condensate as well as its Taylor coefficients. Since we work with staggered quarks, the mass renormalization is 
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FIG. 1: A diagrammatic representation of the terms in the expansion of the chiral condensate, extending the rules in [l3jl . 
The lines represent quark propagators. The dots are insertions of M' (which is 70 in the continuum limit) and the crosses 
are mass insertions. They represent derivatives with respect to fi and m respectively. The contribution of each topology to 
the derivatives in eq. 1131 can be found by simple counting. Every closed loop is a Fermion trace. These diagrams should be 
thought of as being connected to all possible topologies of gluons. 



multiplicative, m R (a) = V Z m (a)m, where the renormalization constant depends on the ultraviolet cutoff, the lattice 
spacing a, but not on T and \i. Throughout this paper we write the bare quark mass as m. The quark condensate is 
also renormalized by Z m . If we are not interested in the absolute value of the condensate, but only on its variation 
with [i at constant T, then one can simply work with either of the ratios 

Cs{m R ,T,n u ,n d ;a) Cs(m R ,T, fj, d ;a) . 

— — Z scheme or — — T scheme . 12 

Cs(m R , 0,0,0; a) CsymR, T, 0, 0; a) 

One expects the T scheme to be undefined for T > T c when m R = 0, since the chiral condensate then vanishes. 
Elsewhere there is a finite multiplicative relation between the schemes which has the interpretation of the T-dependence 
of Cs at fj, = 0. The computations have to be performed at fixed renormalized mass. We perform a formal expansion 
here which will later be renormalized in either of these schemes. 

For the u quark condensate (C u ) — (uu) — (RetrM" 1 ), clearly C« = ci" d ' ) = Cif — 0. The non- vanishing 

derivatives are d u) = -Retr M- l M' u M~ l = C, and C^ u) = Retr (2M~ l M^M' 1 M' U M- X - M^ l M'^M~ l ) = C" . 
For the d quark condensate, the corresponding relations can be found by flipping u and d indices. Since the Taylor 
series expansion is carried out at fj, = 0, in the limit of degenerate masses the quantities C and C" are flavour 
independent. In this limit the Taylor expansion coefficients of the flavoured condensates are flavour independent, and 
we need to evaluate the four quantities — 

C° s = (C), 

Cl = {20 X C +C), 

Cf = (2(0 11 +0 2 )C + 20 1 C'+e")-2(O u +0 2 )(C), 

Cl 1 = 2(O n C + 1 C')-2(O n )(C), (13) 

A direct check shows that Cg vanishes, as does every odd term in the expansion. C|° is the diagonal, and Cg 1 is the 
off-diagonal QRC. 

A diagrammatic representation of the expansion of the chiral condensate is shown in Fig. ^ This builds on a 
previous diagrammatic representation of the expansion of the partition function |13( . Every line denotes a fermion 
propagator, a dot denotes a derivative with respect to /// (an insertion of 70 in the continuum), every cross is a 
derivative with respect to (a mass insertion), and every closed line is a trace over Fermion coordinates. The 
fermions have to be dressed in all possible ways by gluon insertions. Power counting arguments |15| can be easily 
extended to estimate the power of the gauge coupling g associated with the connected parts which appear in eqs. i|T3)) 
when they are evaluated at high temperatures. The largest quark-loop-disconnectcd quantity is {OiC) c which is of 
order g 4 . Both (On) and (0\C') C are of order g 6 , and (Oi\C) c is of order g 8 . These powers may be modified by logg 
terms due to infrared singularities. (C") is of order g°, and, because it is non- vanishing, generically gives the largest 
contribution in the high temperature phase. 
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The Taylor expansions of Cg are 

C s (m, T, M3 ) - + (Cf° - C* 1 ) ^ + • • • (14) 

where the temperature dependence of the coefficients has not been shown explicitly. Note that in the T scheme the 
Taylor series begin with a constant term of exactly unity. In this scheme one can rewrite the series as an expansion 
in \ijT and render all the coefficients dimensionless. Then the remaining Taylor coefficients have the significance of 
fractional changes due to non-zero values of /i. The power counting argument shows that for T ^> T c the coefficient 
of the quadratic term is dominated by the operators (C") c , and the two expansions differ only by terms subleading in 
g. In the low-temperature phase these arguments fail, and the two series could have significantly different behaviour. 
Near and below T c other arguments, outlined in Section llVl also lead us to expect different behaviour for the two 
series. 



C. Chiral Ward identities for meson susceptibilities 

A meson operator is Af(x) = ■if;(x.)Tip(x.), where ip is a quark field operator and T is a Hermitean spin-flavour 
matrix. The conjugate operator is M^(x) = ip^T^jotp — ipTip — ±ipTip, where the plus sign is obtained for the Dirac 
matrices 1, 70 and 7570 and the negative sign for 75, 7* and 75 7;. The meson correlators are 

C(x,y) = (M(x)Mt(y)) = (trrAf- 1 n(y)f(Mt)- 1 n(x)), (15) 

where M is the Dirac operator, and II(x) is the projector to one site x. The final expression on the right was obtained 
using the relation II 2 = II and the fact that II and T commute since they act on different spaces. In Q the correlator 
was parametrized in the form A cosh m(L/ 2 — x) and the derivatives with respect to /i of the parameters A and m were 
related to operator expectation values obtained by taking derivatives of the right hand side of eq. (|15|) . However, since 
CP symmetry is explicitly broken on introduction of chemical potentials, one need no longer have C(x, y) = C(|x— y|). 
This is explicitly visible in the computations reported in |16l Il7| . 

A little tinkering does not help, since it is the transfer matrix formalism that needs to be re-examined at finite /x, 
because a non-zero value of this parameter introduces an explicit breaking of time reversal symmetry. It is much more 
straightforward to deal with the meson susceptibilities 1 

XT 4E C ( X ^) = i<trrM- 1 r(Mt)- 1 ). (16) 

If a single pole dominates this propagator then the corresponding mass, mr cx 1/^/xr- However, this quantity more 
generally measures mesonic fluctuations |lq. 

Chiral Ward identities are consequences of certain operator equalities which follow from chiral symmetries. Each 
formulation of lattice fermions can preserve a different subset of these identities; a basic set for staggered fermions 
can be found in [l9| . The remaining continuum identities are violated at finite lattice spacing, and hence can be 
used to check the approach to the continuum. The prototypical chiral Ward identity is Cs(T,fi) = m% 7r (T, /1), where 
X-n is the pseudo-scalar susceptibility. Furthermore, this chiral Ward identity allows us to define a renormalized x-k 
using the simultaneous renormalization of the quark mass and the condensate. Taylor expansions of both sides of this 
identity can then be equated term by term. Another conclusion is that the Taylor expansion of Xtt must be even in 
/i. Explicit construction of the first Taylor coefficient shows that it vanishes, in conformity with this proof. 

A second chiral Ward identity [2(j relates the isovector scalar susceptibility to the mass derivative of the condensate: 
dCs/dm = — Xe- Again this can be used to renormalize Xe an d relate two Taylor series. We emphasize that this 
involves the isovector scalar, i. e., the scalar meson susceptibility implicated here is obtained out of a connected quark 
loop. 

We show later that in the continuum limit of the high temperature phase the second derivatives of the condensate 
with respect to the chemical potential vanish. As a result, the QNS arc insensitive to variations of the quark mass, 
and the pion susceptibilities are insensitive to chemical potential. Also, in this phase, symmetry arguments show that 
Xtt = Xe H3 > implying that Cs vanishes linearly with to. As a result, we can deduce that the scalar susceptibility 
is also independent of the isoscalar chemical potential, at least to quadratic order. Below T c this chain of logic does 
not hold. The scalar susceptibility is interesting because of speculation about the massless modes at the critical end 
point |22(. We will discuss it at greater length elsewhere. 
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D. Maxwell relations 

A Maxwell relation is the equality of two distinct physical interpretations of a mixed derivative obtained by inter- 
changing the order of the derivatives. From the Taylor expansion of the chiral condensate in eq. (|14|l we can find 
Maxwell relations with the change of quark number susceptibilities with the quark mass. The leading order relation 

dfi dm 

was first noted in p3 |. It is trivially true at fi = since the first derivative on the left vanishes, as does n for all quark 
masses. 

The second derivatives give two non-trivial Maxwell relations and consequent relations using the chiral Ward 
identities discussed earlier, 



C 



2(1 



— m- — 



s dfi u dfj,d d^ud^id dm 

There are higher order relations between the mass derivative of the non-linear susceptibilities and the higher Taylor 
coefficients of the condensate. This class of Maxwell relations can be checked by explicit differentiation. Preliminary 
results have been reported in and a mixed version of this has been used later to correct for the mass dependence 
of a Taylor expansion of the pressure [25J. Here, as a byproduct of our computation of the renormalized chiral 
condensate, we shall give the continuum limit of the derivative of the susceptibilit y. 
The relative rates of strange and light quark production in heavy-ion collisions |2y| is 

A s = ^i, (19) 

Xuu 

with obvious extensions to the production rates of heavier quarks. Since it is hard to perform lattice computations at 
realistic values of light quark masses due to constraints of computer time, one can lighten the computational burden 
by using a Taylor series for the mass dependence of the susceptibilities utilizing the above Maxwell relation. Then 
one can compute A s at some reasonably light quark mass, corresponding to, say, the pion mass being two to three 
times heavier than in the real world, and extrapolate to the physical quark mass values using the Maxwell relation. 
One can also investigate the stability against small variations in the heavy quark mass by similar means. 



III. RESULTS 



A. Implementing the renormalization scheme 

In eq. (|12[1 we have defined a renormalization scheme in which the multiplicative renormalization of the chiral 
condensate at temperature T and chemical potential fi is absorbed by the bare condensate at T = or /i = 
evaluated at the same lattice spacing a. In order to implement this scheme this compensation must be performed at a 
constant renormalized quark mass. We complete the renormalization prescription by choosing to trade the pion mass 
for the renormalized quark mass. The practical question is then the choice of the bare quark mass which corresponds 
to a fixed pion mass. 

The lattice spacing has been connected to the temperature scale using the QCD /3-function at full two-loop order 
[2?j |. This scale determination was used to tune the bare couplings in the Wilson action required to simulate the 
quenched theory at a fixed physical temperature while changing the lattice spacing [2^| . Here we used the same bare 
couplings to extract the pion mass using T = simulations. 

The series of bare couplings used correspond to fixed T = 1.5T C on lattices with Nt — 4, 8, 10 and 12 and to fixed 
T = 0.75T C for Nt — 8, 10 and 12. The lattice sizes, L, were selected so that Lm, > 6 in order to keep finite volume 
effects under control. At the two coarsest lattice spacings a previous determination of meson masses [29| could be 
used to estimate m v and select an appropriate value of L. Then one could successively use the already determined 
masses to estimate the size which would be required at finer lattice spacing. Quenched configurations were generated 
using the Wilson action. At each lattice spacing the first 500-1000 sweeps were discarded for thermalisation and a 
configuration was kept after every subsequent 50 sweeps. 
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aT c 




size 


am q 


2 2 

a 


x 2 


Range 


am p 


2 

X 


Range 




0.1667 


5.8941 


16 3 x 32 


0.01667 


0.119(1) 


1.99 


7-10 


0.66(2) 


1.43 


5-17 


0.52(2) 








0.02 


0.141(2) 


2.82 


7-10 


0.68(2) 


1.72 


5-17 


0.55(2) 








0.03 


0.205(2) 


1.39 


8-10 


0.74(2) 


2.54 


5-17 


0.61(2) 






30 4 


0.01667 


0.123(3) 


4.41 


6-11 


0.71(4) 


10.73 


3-16 


0.49(3) 








0.02 


0.142(2) 


5.59 


6-10 


0.72(3) 


11.08 


3-16 


0.52(2) 








0.03 


0.214(4) 


4.10 


6-11 


0.77(3) 


11.78 


3-16 


0.60(2) 


0.1333 


6.0150 


20 3 x 40 


0.01333 


0.080(1) 


1.66 


7-13 


0.56(2) 


4.95 


5-21 


0.50(2) 


0.1111 


6.14 


24 3 x 48 


0.01111 


0.0515(9) 


1.64 


9-13 


0.42(1) 


2.62 


7-25 


0.54(1) 


0.0833 


6.3384 


40 4 


0.0083 


0.035(1) 


6.58 


9-17 


0.333(9) 


2.55 


7-21 


0.56(2) 








0.01 


0.040(1) 


5.06 


10-17 


0.341(8) 


2.64 


7-21 


0.59(2) 


0.0667 


6.525 


42 4 


0.0067 


0.022(1) 


2.45 


13-19 


0.25(1) 


7.35 


8-22 


0.60(3) 








0.005 


0.018(1) 


2.85 


13-19 


0.25(1) 


6.67 


8-22 


0.54(3) 


0.0556 


6.65 


50 4 


0.0056 


0.0161(3) 


8.13 


13-23 


0.21(1) 


15.0 


10-25 


0.60(3) 



TABLE I: Meson masses in quenched QCD with Wilson action and staggered quarks. The \ 2 f° r pions is obtained by fitting 
local masses to a constant. The rho meson mass is obtained with a four parameter fit to the vector correlator. Details are 
discussed in the text. 
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FIG. 2: The dimensionless ratio m^(T = 0)/mT for fixed bare quark mass m/T c = 0.1 at a succession of lattice spacings 
corresponding to fixed T = 1.5T C at varying N t shown as a function of l/N? — a 2 T 2 . The constancy of the ratio shows that 
this procedure correctly holds the renormalized quark mass constant and thereby implements the renormalization scheme of 



eq. 



The stored configurations were used to determine the chiral condensate, through a noisy evaluation of the trace 
of the inverse of the fermion matrix, and meson correlators. The stopping criterion on the staggered fermion matrix 
inversion was decreased by factors of 10 until the extracted value of m ff was seen to stabilize. It was observed that 
while decreasing lattice spacing keeping fixed the stopping criterion had to be decreased significantly. We attribute 
this to the fact that as the lattice spacing decreases at zero temperature the eigenvalues near a rapidly decreasing 
renormalized quark mass take more conjugate gradient steps to resolve. 

Our extraction of pion masses was based on effective (local) mass extraction, m(t), and cross checked by fits. As 
has been known for a long time, this is contaminated at small distances by higher lying states and at large distances 
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0.005 




FIG. 3: The diagonal QRC in the T scheme, at fixed T/T c = 1.5 holding the renormalized quark mass fixed by fixing the bare 
mass to be m/T c = 0.1. The full lines show the 1-cr band on the continuum extrapolation and the dashed lines show the 3-cr 
confidence band. 
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FIG. 4: The diagonal QRC in the T scheme, at fixed T/T c = 0.75 holding the renormalized quark mass fixed by fixing the bare 
mass to be m/T c — 0.1. The lines show the 1-cr band on the continuum extrapolation. 



by periodic lattice artifacts. When a plateau can be seen, that value is taken as an estimate of the ground state mass. 
We converted this into an automatic task by varying i min and i max and for each fixed value of this parameter checking 
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whether the data on m(t) could be fitted to a constant. We accepted the largest interval over which this could be 
done and took the corresponding constant as an estimate of the mass (see Table Q] For the vector meson the local 
mass is harder to stabilize, and we used four parameter fits to extract the mass, varying the fit interval to look for 
stability. 

In simulations at the coarser lattice spacings we performed a series of computations at larger bare quark mass. These 
are useful in order to check that the lattice volumes are under control. The near constancy of the ratio a 2 m^ r /(am q ) 
at fixed a on the coarser lattices indicates that the lattice volumes are sufficiently large that one can extrapolate to 
the correct continuum physics. 

In simulations at a sequence of lattice spacings, a, corresponding to a fixed temperature and varying N t (a = l/TN t ), 
with a sequence of bare quark masses m/T c = 0.1, we found that the dimcnsionless ratio m^/rriT is fixed. This is 
shown in Figure [21 The constancy of this ratio shows that this sequence of bare quark masses corresponds to a fixed 
renormalized quark mass, and gives an easy method to tune the bare quark mass in order to fix the renormalized 
quark mass. At this quark mass, we found that m w /m p 0.6, indicating that m p could also have been chosen to fix 
the renormalized quark mass. For the parameters chosen here, the fact that the renormalized quark mass is rather 
heavy is seen both in the fact that the ratio m^/mp is larger than the physical value and also in that m^/Tc « 2.3. We 
recall that with two flavours of dynamical staggered quarks choosing m/T c — 0.1 corresponds to m w /m p = 0.36±0.01, 
and that on reducing the bare quark mass to m/T c = 0.03 the ratio m^/Tc became physical |30|. 

B. Continuum limit in quenched QCD 

1. T>T C 

We used stored gauge configurations from the study in [2^ for our measurements. These were obtained on Nt — 
4,8, 10, 12 and 14 lattices for temperatures T = 1.5T C ,2T C and 3T C respectively. The main concern in choosing the 
spatial size N s is to avoid gross distortions of the infinite volume result. For this it is sufficient to have N s /N t > 
max(2,T/T c ). For Nt = 4 we have used N s = 12. For all the other Nt the spatial size was chosen to be N s = 2N t + 2 
for T/T c — 1.5 and 2, and N„ = 3N t + 2 for T/T c — 3. For details of the generation of the configurations and their 
statistical properties, we refer the reader to |28j . In this work we recomputed, as a check, the operators 0±, On and 
02 which were studied in [2^| and supplemented these by computations of C, C and C" needed for the computation 
of the Taylor expansion of the chiral condensate. The traces were evaluated by the noisy method (see [28| for details). 
For the conjugate gradient algorithm to obtain the matrix inverses, we used a stopping criterion that the norm of the 
residual should be less than 10~ 3 y/ NtNf . 

The vanishing of the expectation values (Oi) and (C) are checks of the accuracy of the computation. On each 
lattice and at every temperature, we found that the expectation values were zero to within a part in 10 6 with the 
above stopping criterion. In the high temperature phase of QCD, where we work, Xud is small, and therefore one 
expects a noisy and almost vanishing signal for its derivative with respect to the quark mass. In agreement with this 
expectation, we find that Cg 1 is too small and noisy to be measured. We return to this in Section HVI 

The diagonal QRC in the T scheme is shown in Figure |21 It is clear that there is significant /^-dependence on 
coarser lattices; in fact for N t = 4 we found that T 2 Cf{T)/C s {T) ~ -0.096 at T/T c = 1.5. However, this Taylor 
coefficient vanishes at the 99% confidence limit on extrapolation to the continuum. The conversion from T scheme 
to Z scheme involves a finite multiplicative factor Cs(mn^ T)/Cs(rriR, 0) which also has an interpretation as the 
temperature dependence of the condensate. At T = 1.5T C we find on extrapolation to the continuum that this factor 
is 0.26 ± 0.01 for ma such that m n /m p ~ 0.6. With decreasing quark mass this factor is expected to decrease. 

A vanishing continuum limit of the QRC is also found for m/T c = 0.1 at the two other temperatures that we have 
studied. The identities in eq. (|18|) then imply that Xuu is insensitive to changes in the quark mass and x-n is insensitive 
to fi in this range of temperatures. Since Xuu agrees with a perturbative evaluation for T > 1.5T C , its insensitivity to 
m can be understood from the fact that the effective infrared cutoff is given by the Matsubara frequency 7rT and not 
by m when m/T c < ttT/T c . 

2. T <T C 

We also made a series of simulations at fixed T/T c = 0.75 for N t = 6, 8, 10 and 12 in quenched QCD with 
the Wilson action. Three hits of pseudo-heat-bath followed by one Metropolis step made up one composite sweep. 
The first 500 sweeps were discarded for thermalization and subsequently one configuration was saved once every 500 
sweeps. Between 35 and 45 such stored configurations were used on each lattice for this work. Measurements were 
made with the bare quark mass m/T c = 0.1 and with 50 noise vectors for evaluation of quark traces. In all cases 
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FIG. 5: The diagonal QRC as a function of the temperature in Nf = 2 QCD within the Z scheme, for a — 1/4T holding the 
renormalized quark mass fixed by fixing the bare mass to be m/T c = 0.1. Pluses denote results obtained on 4 x 8 3 lattices, 
pentagons on 4 x 12 3 and circles on 4 x 16 3 . The boxes represent the quenched QCD value. 




FIG. 6: The off-diagonal QRC as a function of the temperature in Nf = 2 QCD within the Z scheme, for a = 1/4T holding the 
renormalized quark mass fixed by fixing the bare mass to be m/T c — 0.1. Pentagons represent results for 4 x 12 3 , and circles 
for 4 x 16 3 . 



we used spatial sizes N s = 2N t . From the T = runs described earlier, it can be seen that the size of the spatial 
box corresponds to over 5.5 pion Compton wavelengths (at T — 0) and hence finite size effects are under reasonable 
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FIG. 7: The continuum extrapolation of the diagonal QRC renormalized in the Z scheme in Nf = 2 QCD with bare quark 
mass of 0.1T C . The circles denote data for T = 0.9T C , and the 3-cr band for the continuum extrapolation is shown. The boxes 
represent data at T c , and the 1-a (nearer arrows) and 3-cr (further arrows) limits of the extrapolation shown are discussed in 
the text. 



control. 

We found that the continuum limit of both (O2) and (On) were consistent with zero within reasonably small errors. 
On all but the coarsest lattice, N t — 6, the first expectation value was significantly smaller than the second. As a 
consequence Xud and Xuu were equal within errors even at finite lattice spacing. As shown in Figure ^ Cl also 
vanishes within errors in the continuum limit, showing that the vanishing of Xuu is not an accident of the choice of 
the quark mass. The off diagonal QRC was consistent with zero at each N t . 



C. Dynamical staggered quarks 



We report the results of a series of computations with two flavours of dynamical quarks as the temperature varies 
between 0.75T C and 2T C . The simulations were performed using the R- algorithm with two flavours of quarks with bare 
mass of 0.1T C , corresponding to m^/mp = 0.36 ± 0.01 30] . Three lattice sizes were used in most of the simulations, 
namely 4 x 8 3 , 4 x 12 3 and 4 x 16 3 . The critical coupling 8 C for Nt = 4 was determined in |3l| and the setting of the 
temperature scale was performed by methods given in |3J| ■ The critical coupling for this quark mass was seen in the 
peaks of Xe, the Wilson line susceptibility, and even in autocorrelation times |30l l3l| . The critical coupling is also 
consistent with those found for neighbouring values of quark mass in 

Lattice computations with dynamical quarks are still too costly for the continuum limit to be taken easily. However, 
in order to have some idea of the continuum limit we have also performed two further simulations with 6 x 12 3 and 
8 x 16 3 lattices at T = 0.9T C . In addition, we have results at T c on 6 x 12 3 lattices. The critical couplings at these Nt 
were determined in [33^ . It is impossible to collect a sufficient number of statistically independent configurations on 
N t = 8 lattices near T c with current computing resources due to a rapid increase in autocorrelation times, r. In all 
our simulations autocorrelations of local observables were monitored and the total run length was at least 50 times 
the longest, r max , when measured self consistently over the history. The measurements used between 50 and 250 
configurations, each separated by r max . We work in the Z scheme since we are interested in the T dependence of the 
derivatives. 

Consistent with our previous observations, Cg was seen to vanish with small errors across the full range of tem- 
peratures considered. The diagonal QRC was negative and differed significantly from zero, as shown in Figure [3] For 
comparison, results from the quenched theory at the same bare quark mass are also shown. The T-dependence is 
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FIG. 8: The panels on top show the expected phase diagrams in the T-/j,o and T-/X3 |j| planes for m,R > 0. The panels below 
show the expectations in the chiral limit, tur = Hlf . Thick lines denote second order transitions, the tri-critical point is 
marked with a fancy circle and the dotted line shows a part of the critical line which could be a model artifact |3fij . Further 
interesting phases at large fi are not shown. 



seen to be a little different; this could reflect either the difference in the value of m^/nip or be a quenching artifact. 
The off-diagonal QRC (see Figure [SJl is consistent with zero within small errors above T c . Near and below T c larger 
fluctuations are seen. Although still consistent with zero at the 3-cr level, the averages increase by several orders of 
magnitude compared to T > T c , and seem to be comparable in magnitude to the diagonal QRC. 

The results of our preliminary investigation of the continuum limit with dynamical quarks are shown in Figure \7\ 
At 0.9T C the 3-cr error band on the extrapolation of the diagonal QRC is consistent with zero. At T c the continuum 
extrapolation perforce had to be performed with data from only two lattice spacings. These completely fix the 
parameters of a linear extrapolation to the continuum. Error estimates were obtained by then taking the slope as this 
given parameter and estimating the error on the continuum value in the usual way. The resulting l-a and 3-cr error 
limits on the extrapolation are shown in Figure Another estimate is to take the extrapolations of the upper end of 
one error bar and the lower end of the other. This lies within the 3-cr band given. With all these caveats, our data 
currently points to a non-zero continuum extrapolation of renormalized Cg° at T c . 



IV. DISCUSSION 



The chiral condensate is a derivative of the free energy density with respect to quark mass (eq. 1111) . and the 
latter requires renormalization. As a result, so does the condensate and its Taylor coefficients. We have defined two 
renormalization schemes for staggered quarks (eq. in order to set up a well defined Taylor series expansion for the 
chiral condensate as a function of the chemical potential /j, (eas.H3land ll4|) . Diagrammatic rules for writing down the 
operators defining the Taylor coefficients have been given (Figure^. The first order Taylor coefficient, C s vanishes by 
symmetry. The two QRCs are related through Maxwell relations and chiral Ward identities to all other second order 
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Taylor coefficients in the chiral sector feas. 1180 . The diagonal QRC, C|°, has small errors and is easily measurable 
at the renormalizcd quark mass we used. We found that at finite lattice cutoff it is negative. Away from T c , both in 
quenched and Nf — 2 QCD, the continuum limit of this quantity vanishes. We have presented evidence (Figure EJ) 
that near T c in the Nf — 2 theory it does not vanish. 

This negative peak in C|° (Figures |5] and [JJ corresponds to a significant shift of the condensate, and hence Xki 
with chemical potential. The implications for strangeness production have been pointed out in Section III Dl Here we 
discuss how measurements of the diagonal and off-diagonal QRCs, C|° and Cg 1 , test current conjectures about the 
phase diagram of QCD at small \i based on the perturbative expansion or model computations |a, |9l llCt lll| . 

The location of T c is often identified through a peak in the scalar susceptibility \e (which is the same as the mass 
derivative of the chiral condensate through the Patel Ward identity HJ})- In the chiral limit this is an unique definition, 
since the direction of spontaneous symmetry breaking at low T makes this a massive particle, which becomes massless 
precisely at T c . However, at finite quark mass the symmetry is explicitly broken, Yjl is finite throughout the low T 
phase and there is a cross over instead of a critical point along the T axis at T c |36j, as shown in the top panel of 
Figure|Hl(see however, 37]). The cross over temperature T c is usually not uniquely denned, in the following sense. At 
a phase transition different physical quantities give exactly the same estimate of T c . This is no longer true at a cross 
over. Thus the slight shift in the peak of Cg (Figure |SJ) compared to T c previously defined from peaks in Xe and the 
Wilson line susceptibilities |3(J, |31j lends some support to the scenario of a cross over. 

Furthermore, the observed negative peak in the quadratic response of the condensate shows that x™ decreases 
at finite /i, and hence pions cannot be the soft mode as one approaches the critical end point. To the best of our 
knowledge, this is the first such indication from lattice QCD. Unless C|° changes sign in the continuum, this argument 
continues to hold. We have presented evidence that it does, because the diagonal QRC does not turn positive in the 
vicinity of T c (Figure EJ. Thus support for the standard picture of the phase diagram for vtlr > is provided by C s 
in two different ways — by indicating that at fi = there could be a cross over, and by showing that pion fluctuations 
decrease as the baryon density increases. 

It is interesting to see what else the QRCs can say about the phase diagram. For tur > 0, as T decreases, the 
Taylor series for Cs will eventually hit the critical end point in an extrapolation in /iq |8J|, and the critical line in /X3 
|^.l34| (see upper panels in Figure[SJ). Therefore, at lower T there must be differences in the two extrapolations. This 
can only happen if the off-diagonal QRC becomes significant. We have shown preliminary evidence for this in Figure 
|SJ This should become easier to measure with decreasing vtlr since m T as well as the chemical potential at the critical 
end point are expected to decrease. 

In the chiral limit the phase diagram in the T-^o plane is expected Q to contain a critical line meeting the T axis 
at a critical point T c . It is also expected that this critical line ends in a tri-critical point where it meets a line of first 
order phase transitions emerging from the hq axis (see lower panels in Figure The phase diagram in the T-fx^ 
plane has been explored in lid I 111 . A second critical line is expected to emerge from the T axis, curving down to 
eventually touch the 113 axis |35l |. Since there is a unique transition along the T-axis for /io = A*3 = 0, one expects 
both critical lines to start from T c . 

A Taylor expansion cannot be started at the critical point T c , but it is possible to make such a series expansion 
immediately above and below, and thus track the curvature of the two critical lines. If the off-diagonal QRC vanishes, 
then the physics along these two lines is the same. However, one has massless pions, whereas the other does not. 
Also, at lower temperature, one series should give evidence of a tri-critical point on extrapolation, whereas the other 
should not. Thus, if the phase diagrams drawn above are accurate, then they should be reflected by a non- vanishing 
value of the off-diagonal QRC. Further, measurements of Cf° ± C5 1 would give information on the soft modes in these 
two directions. 

Thus, there are several extensions of this study which would be interesting to perform in the future. First, the chiral 
extrapolation of Cl and Cg 1 obtained at fixed N t by varying ttir would indicate whether the tri-critical point indeed 
appears at tur = 0. If the standard picture does hold, then the chiral extrapolation of C|° + Cg 1 and C|° — Cg 1 at 
T just below T c would give information on the soft modes along the two critical lines. 

It is a pleasure to acknowledge discussions with Rajiv Gavai, Edwin Laermann and Maria-Paola Lombardo. Part 
of the computations were carried out on the Indian Lattice Gauge Theory Initiative's CRAY XI at the Tata Institute 
of Fundamental Research. 
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